#############################
#############################
###Spatial Regression Code###
#############################
#############################
library(sp)
library(spdep)
library(sf)
library(plm)
library(splm)
library(spatialreg)
####Subset data without NAs
subset.dat3 <- na.omit(full.1996.2019[,c("acled_battle_state","acled_battle_rebel","acled_battle_polmil","acled_battle_idmil","latitude","longitude","log.nl","logged.pop","lagacled_battle_state","lagacled_battle_rebel","lagacled_battle_milall","lagacled_battle_polmil","lagacled_battle_idmil","t","spei","p_anom","outbreak","ym","t_anom","gid","month", "year", "COWCODE")])  
####Subset years of analysis 
subset.dat3 <- subset(subset.dat3, year>=2000&year<=2018)
###Create a balanced panel using plm
afro.plm<- pdata.frame(subset.dat3,index=c("gid","ym"))
afro.plm.gid<- make.pbalanced(afro.plm,balance.type = "shared.individuals")
afro.plm.balanced<- make.pbalanced(afro.plm.gid,balance.type = "shared.times")

##Check number of cells
length(unique(afro.plm.balanced$gid))

###Import polygon data and create relevant measures
polygons.prio<- sf::st_read("priogrid_cellshp/priogrid_cell.shp")
polygons.prio<- polygons.prio %>% dplyr::select(gid)
polygons.prio.slice<- subset(polygons.prio, gid %in% afro.plm.balanced$gid)
#rownames(polygons.prio.slice)<- 1:4064

##Create a GID indicator corresponding to the polygons
gid<- afro.plm.balanced
gid$gid <- as.numeric(paste(as.factor(gid$gid)))
gid$one <- 1
gid <- summaryBy(one~gid, data=gid, keep.names = T, FUN=c("mean"))
gid$indiv <- (1:length(unique(gid$gid)))
gid$one <- NULL
summary(gid)

##Create spatial weights
slice.sp<- as(polygons.prio.slice,"Spatial")
slice.nb<- poly2nb(slice.sp,queen=T)
# Row standardized
slice.listw<- nb2listw(slice.nb,style = "W",zero.policy = T)

#If needed, index variables by the first two columns
#afro.plm.balanced<- merge(afro.plm.balanced,gid,by="gid")

#Reformat dataset as a plm object if needed
#afro.plm.balanced<- pdata.frame(afro.plm.balanced,index=c("gid","ym"))
#Create spatial outbreak lag
afro.plm.balanced$s_outbreak<- slag(afro.plm.balanced$outbreak,listw = slice.listw)
#memory.limit(size=1000000)


###############################
###Pooled Spatial Lag Models###
###############################
#####State conflict
cellFE.sldv<- spml(acled_battle_state ~ 
                     log.nl + logged.pop +
                     lagacled_battle_state + t  +
                     p_anom+spei+outbreak+t_anom+s_outbreak+
                     factor(month),
                   data=afro.plm.balanced,listw=slice.listw,na.action=na.omit,
                   index=c("gid","ym"),spatial.error = "none",
                   lag=T,model="pooling")
summary(cellFE.sldv)

#####Rebel conflict
cellFE.sldv2<- spml(acled_battle_rebel ~ 
                      log.nl + logged.pop +
                      lagacled_battle_rebel + t  +
                      p_anom+spei+outbreak+t_anom+s_outbreak+
                      factor(month),
                    data=afro.plm.balanced,listw=slice.listw,na.action=na.omit,
                    index=c("gid","ym"),spatial.error = "none",
                    lag=T,model="pooling")
summary(cellFE.sldv2)

####Political militias
cellFE.sldv3<- spml(acled_battle_polmil ~ 
                      log.nl + logged.pop +
                      lagacled_battle_polmil + t  +
                      p_anom+spei+outbreak+t_anom+s_outbreak+
                      factor(month),
                    data=afro.plm.balanced,listw=slice.listw,na.action=na.omit,
                    index=c("gid","ym"),spatial.error = "none",
                    lag=T,model="pooling")
summary(cellFE.sldv3)

####Identity militas
cellFE.sldv4<- spml(acled_battle_idmil ~ 
                      log.nl + logged.pop +
                      lagacled_battle_idmil + t  +
                      p_anom+spei+outbreak+t_anom+s_outbreak+
                      factor(month),
                    data=afro.plm.balanced,listw=slice.listw,na.action=na.omit,
                    index=c("gid","ym"),spatial.error = "none",
                    lag=T,model="pooling")
summary(cellFE.sldv4)


##############################
###Random Effect Lag Models###
##############################
#####State conflict
cellFE.sldv.re<- spml(acled_battle_state ~ 
                        log.nl + logged.pop +
                        lagacled_battle_state + t  +
                        p_anom+spei+outbreak+t_anom+s_outbreak+
                        factor(month),
                      data=afro.plm.balanced,listw=slice.listw,na.action=na.omit,
                      index=c("indiv","ym"),spatial.error = "none",
                      lag=T,model="random", effect="individual")
summary(cellFE.sldv.re)

#####Rebel conflict
cellFE.sldv2.re<- spml(acled_battle_rebel ~ 
                         log.nl + logged.pop +
                         lagacled_battle_rebel + t  +
                         p_anom+spei+outbreak+t_anom+s_outbreak+
                         factor(month),
                       data=afro.plm.balanced,listw=slice.listw,na.action=na.omit,
                       index=c("indiv","ym"),spatial.error = "none",
                       lag=T,model="random", effect="individual")
summary(cellFE.sldv2.re)

####Political militias
cellFE.sldv3.re<- spml(acled_battle_polmil ~ 
                         log.nl + logged.pop +
                         lagacled_battle_polmil + t  +
                         p_anom+spei+outbreak+t_anom+s_outbreak+
                         factor(month),
                       data=afro.plm.balanced,listw=slice.listw,na.action=na.omit,
                       index=c("indiv","ym"),spatial.error = "none",
                       lag=T,model="random", effect="individual")
summary(cellFE.sldv3.re)

####Identity militas
cellFE.sldv4.re<- spml(acled_battle_idmil ~ 
                         log.nl + logged.pop +
                         lagacled_battle_idmil + t  +
                         p_anom+spei+outbreak+t_anom+s_outbreak+
                         factor(month),
                       data=afro.plm.balanced,listw=slice.listw,na.action=na.omit,
                       index=c("indiv","ym"),spatial.error = "none",
                       lag=T,model="random", effect="individual")
summary(cellFE.sldv4.re)
